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The dynamics of macroscopically homogeneous sheared suspensions of neutrally buoyant, 
non-Brownian spheres is investigated in the limit of vanishingly small Reynolds numbers 
using Stokesian dynamics. We show that the complex dynamics of sheared suspensions 
can be characterized as a chaotic motion in phase space and determine the dependence of 
the largest Lyapunov exponent on the volume fraction <f>. We also offer evidence that the 
chaotic motion is responsible for the loss of memory in the evolution of the system and 
demonstrate this loss of correlation in phase space. The loss of memory at the microscopic 
level of individual particles is also shown in terms of the autocorrelation functions for the 
two transverse velocity components. Moreover, a negative correlation in the transverse 
particle velocities is seen to exist at the lower concentrations, an effect which we explain 
on the basis of the dynamics of two isolated spheres undergoing simple shear. In addition, 
we calculate the probability distribution function of the velocity fluctuations and observe, 
with increasing 0, a transition from exponential to Gaussian distributions. 

The simulations include a non-hydrodynamic repulsive interaction between the spheres 
which qualitatively models the effects of surface roughness and other irreversible effects, 
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such as residual Brownian displacements, that become particularly important when- 
ever pairs of spheres are nearly touching. We investigate the effects of such a non- 
hydrodynamic interparticle force on the scaling of the particle tracer diffusion coefficient 
D for very dilute suspensions, and show that, when this force is very short-ranged, D be- 
comes proportional to <j) 2 as <j> — > 0. In contrast, when the range of the non-hydrodynamic 
interaction is increased, we observe a crossover in the dependence of D on </>, from tfi 2 to 
4> as 4> — » 0. 



1. Introduction 

The phenomenon of shear-induced particle diffusion in suspensions at vanishingly small 



Reynolds number has been studied extensively since the work of Eckstein et al. (1977), 
where it was first suggested that a shear flow causes the particles to execute random 
migrations across the streamlines of the ambient flow producing an effect akin to disper- 
sion. This hypothesis of a self-diffusive motion arising from purely viscous hydrodynamic 
interactions between particles has already received considerably experimental support 



(Eckstein et al. 1977; Leighton fc Acrivos 1987; Breedveld et al. 199S, 2001 fljd). More- 



over, as pointed out by Breedveld et al. (2001 fc| ) , the origin of this diffusive behaviour 



is clearly different from the more familiar Brownian diffusion in colloidal suspensions 
caused by thermal fluctuations, as well as turbulent diffusion driven by inertial effects, in 
that this shear-induced diffusion is due only to the hydrodynamic interactions between 
the particles comprising the suspension. Since, in principle, these interactions consti- 
tute a deterministic process, the question arises as to whether and how it can also be 
viewed as a diffusion process. To address this question, it is generally assumed that the 
arrangement of neighboring suspended spheres leads to a time-dependent random event 
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(Eckstein et al. 1977) in that, upon collision with their closest neighbors, the spheres 
will suffer many successive random displacements ultimately leading to a random walk 



(Leighton & Acrivos 1987; Zarraga & Leighton 1999). Underlying this description is the 
fundamental assumption that collisions between spheres eventually become statistically 
independent. Although many theoretical studies have made this strong assumption in or- 
der to calculate the diffusivity of very dilute sheared suspensions from the displacement 
produced by a single collision and then averaging over all initial configurations of the 



colliding spheres ( |Acrivos et al. 1992| ; |Wang et al. 1996| , |1998| ; |da Cunha fc Hinch 1996| ), 
the complex dynamics of suspensions undergoing shear, leading to the loss of correlations 
in the particle motions, has not been investigated thus far in much detail. 

It is the purpose of the first part of this work to pursue the recent suggestion made 



by Marchioro & Acrivos (2001) that the chaotic evolution of sheared suspensions is 
responsible for the loss of memory referred to above, and to demonstrate, via numerical 
simulations, that the evolution of the system in phase space is indeed chaotic. Recall that 
chaotic dynamics is characterized by the sensitivity of the system to initial conditions, 
as evidenced by the exponential growth of the separation distance in phase space of 
two initially neighboring trajectories, and that a standard measure of this sensitivity is 
the largest Lyapunov exponent (LLE), which gives the average rate of this separation 
distance ( Baker fc Gollub 1990| , p. 85). We shall therefore investigate the LLE as a 
function of the volume fraction <fi of the suspensions. We shall also show that the system 
loses the memory of its initial state and that its evolution is asymptotically diffusive. 
The loss of memory at the level of a single sphere will also be discussed in terms of 
the autocorrelation functions of the two transverse velocity components, and we will 
show that, as the concentration is increased and collisions between spheres become more 



4 G. Drazer, J. Koplik, B. Khusid and A. Acrivos 

frequent, the time scale over which the particle transverse velocities remain correlated is 

shortened. 

Before proceeding, it is instructive to consider the motion of a tracer sphere subject 
to purely hydrodynamic interactions with the other spheres in the suspension. As is 
well-known, and as a direct consequence of the linearity of the equations of motion at 
zero Reynolds number, in any encounter between two perfectly smooth spheres neither 
sphere experiences a net lateral displacement, although both may suffer large transient 



displacements from their original streamlines ( Leal 1992 , p. 257). Therefore, for the 
tracer to experience a net displacement leading to diffusive motion it is necessary that 
it interact with at least two other spheres. Since the rate of simultaneous interactions of 
a tracer sphere with two other spheres is proportional to jej) 2 as — > 0, where 7 is the 
shear rate, and since the magnitude of each displacement is proportional to the particle 
radius a, the self-diffusion coefficient should be proportional to jcf) 2 a 2 , in the limit of very 



dilute suspensions (Lcighton & Acrivos 1987). The experimental results by Leighton & 



Acrivos (1987) observe this scaling for volume fractions 0.05 < <p < 0.40. However, in any 
real experiment, suspended particles are not perfectly spherical and, as the separation 
between two colliding particles can be very small (less than 10~ 4 of a particle radius ( |da| 



Cunha & Hinch 1996)), even very small asperities might play a significant role during 



encounters between a pair of spheres. In experiments performed by Rampall et al. (1997 ) 
a roughness of order 10 -3 particle radii was found, and the effect of the asperities on the 
trajectories of nearly touching spheres during a collision was determined. In this situation, 
the interaction of the tracer particle with another sphere will lead to a net displacement 
of the tracer from its original streamline, and therefore, since the rate of interactions 
with another sphere is proportional to 70, the diffusion coefficient is expected to scale as 
70a 2 as 4> — > 0. This linear dependence of the diffusion coefficient on was observed in 
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experiments performed by Phan & Lcighton (1999) and by Zarraga & Leighton (1999) 



On the other hand, thus far, the numerical simulations of sheared suspensions have not 
been extended to low enough values of <f> where one or the other of these two regimes for 
the diffusivity would be expected to apply. 

In this work, we shall therefore investigate the scaling of the diffusion coefficient by 
performing simulations down to values of 4> as low as 0.03. The effect of surface roughness 
and other possible non-hydrodynamic forces, such as residual Brownian forces, will be 
qualitatively modelled by introducing a short-ranged repulsive force between the spheres. 



This is similar to the approach taken by da Cunha & Hinch (199£) and by Zarraga 



& Leighton (1999), where the surface roughness was modelled as a normal force that 
prevents the particles from coming closer than a certain fixed distance. As discussed 
above, the contribution to the diffusive motion of the spheres due to binary collisions 
clearly depends on the magnitude and range of the interparticle force, in that for a 
weak force, a 70 2 a 2 regime for the diffusivity would be expected. However, since at low 
enough concentrations, the linear regime should eventually become dominant no matter 
how small the effect of the interparticle force, we shall investigate this transition from a 
quadratic to a linear dependence of the diffusivity on 4> as <p — > 0. As a preliminary, in 



§2.1 we shall discuss the effects of such a non-hydrodynamic interparticle force on the 
microscopic structure of the suspension, and in particular its dependence on the strength 
and the range of this force. 



2. Simulation method: Stokesian dynamics 

We consider suspensions of non-Brownian particles undergoing shear using the method 
of Stokesian Dynamics which was specifically developed for dynamically simulating the 



behaviour of many particles suspended in a fluid (Brady & Bossis 1984). A detailed de- 
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scription of the method is given in a review by Brady fc Bossis (1988| ), hence only a 



brief discussion is presented here. The method accounts for both hydrodynamic and non- 
hydrodynamic forces between particles. Hydrodynamic forces are computed for spherical 
particles undergoing simple shear, characterized by a shear rate 7, in the limit of zero 
Reynolds number. In simulating the behaviour of infinite suspensions, periodic boundary 
conditions in all direction are imposed, using an adapted version of the Lees-Edwards 



boundary condition in the direction of the shear (|Allen fc Tildesley 1987| , p. 242) ( [Brady 



fc Bossis 1985 ). The volume V of the cubic cell containing a fixed number of spheres N 
is related to the volume fraction <f> by <f> = (Ana 3 /3)N/V. Interactions between particles 
more than a cell apart from each other cannot be neglected, due to the long-range char- 
acter of the hydrodynamic forces, and a lattice sum of the interactions, using the Ewald 



method, is implemented ( Brady et al. 198S| ) 



A typical simulation here consists of N — 64 particles sheared over a period of time 
t ~ IOO7 -1 . The results are averaged over N c ~ 100 different initial configurations using 



the random-phase average method proposed by Marchioro & Acrivos (2001), in order 
to avoid spurious time-periodic fluctuations induced by the time-dependent shape of 
the simulation cell. Each initial configuration corresponds to a random distribution of 
non-overlapping spheres in the simulation cell. 

In what follows, we shall express all the variables in dimensionless units, using the 
radius of the spheres a as the characteristic length and 7 -1 as the characteristic time. 

2.1. Non-hydrodynamic interparticle force 

In a suspension of non-Brownian spherical particles undergoing shear at zero Reynolds 
number, the separation between spheres can be very small (less than 10~ 4 of their radius). 
In this situation, the effects of surface roughness or small Brownian displacements cannot 
be neglected and a short-ranged, repulsive force is usually introduced between the spheres 
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Figure 1. First peak of the pair distribution function g(r) for different values of r c , the charac- 
teristic range of the interparticle force. The position of the peak indicates that, as the range of 
the repulsive force increases, so does the minimum separation between particles. The simulations 
are for (f> = 0.10, N = 64, Fo = 1.0, and N c — 100. g(r) is measured after each initial random 
distribution of spheres is sheared for t ~ 50 and steady state was reached. 

to qualitatively model the behaviour of real systems. The introduction of such a force 
has the further numerical advantage of preventing overlaps between the spheres. 

In this work we used the expression for the repulsive interparticle force, already well- 
tested in the context of Stokcsian dynamics, 

F e~ e/r ' c 

= — l _ c -c/r c e »P ' ( 2J ) 

where 6TT/ia 2 jF a p, with \i being the viscosity of the suspending liquid, is the force exerted 
on sphere a by sphere (3, F n is a dimcnsionless coefficient reflecting the magnitude of 
this force, r c is the characteristic range of the force, e is the distance of closest approach 
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Figure 2. Diagram of a pair of spheres oriented with an angle 8 measured from the downstream 
side of the reference sphere. Down and upstream sides of the reference sphere as well as the 
distance R defining particle pairs are shown. 

between the surfaces of the two spheres divided by a, and e a p is the unit vector connecting 
their centres pointing from (3 to a. 

The magnitude and range of the interparticle force affect the microstructure of the 



suspension (Brady & Morris 1997). As mentioned earlier, colliding spheres in a shear 
flow almost touch one another, with the actual minimum separation distance strongly 



dependent on the repulsion between them. Following Bossis & Brady (1984) we show this 
effect in figure [j] in terms of the pair distribution function g(r), defined as the probability 
of finding the centre of a second particle at a distance r — |r| given that there exists a 
sphere with its origin at r = 0. We see that the minimum separation, and therefore the 
first peak in g(r), is strongly affected by the range of the interparticle force in that, as 
r c increases, the minimum separation between neighboring particles also increases. 

It is also known that the presence of a repulsive force breaks the angular symmetry 
of the microstructure and, in particular, that it destroys the fore-aft symmetry of the 



particle trajectories in a simple shear flow (Bossis fc Brady 1984; Dratler fc Schowalter 



199q). Specifically, again following Bossis & Brady (1984), let us consider the angular 



orientation of pairs formed by spheres closer than a certain distance R, and let us define 
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Figure 3. Normalized angular distribution function gii(9) for pairs of particles. The distance 
between the particles is 2 < r < 2.1 (R = 2.1). Different curves correspond to different values of 
the interparticle force range r c . All simulation were performed with N = 64, iV c = 100, Fo = 1.0 
and 4> = 0.10. 

the probability distribution gu(9) as the probability density of finding a pair having a 



given orientation angle 9 (c. f. figure H). ft is clear that, as already shown by Bossis 



fc Brady (1984 ), particle pairs would be expected to spend more time oriented on the 
upstream side, where the repulsive force is balanced by the shear forces pushing the 
two particles together, implying that the probability of finding a pair oriented upstream, 
90° < 8 < 180°, would be higher than in the downstream range, 0° < 6 < 90°. This is 
borne out by figure || which shows that the pair distribution function becomes increasingly 
more asymmetric, favoring the upstream orientation of particle pairs, as the range of the 
interparticle force is increased. 
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In addition to the microscopic effects mentioned above, the interparticle force also 

affects macroscopic measurable quantities, as we shall discuss in the following sections. 



3. Chaotic motion and time-reversibility 

As already mentioned in the introduction, the question as to whether and how diffusive- 
like transport arises in a suspension of non-Brownian particles undergoing shear has been 



investigated since the original work by Eckstein et aL (1977 ). Experimental evidence 
strongly suggests that even with vanishingly small inertia effects (zero Reynolds num- 
ber) and negligible Brownian and non-hydrodynamic forces, sheared suspensions exhibit 
diffusive behaviour ( Lcighton fc Acrivos 1987 ; Brecdveld et aL 199S , 20016 a|). 

In an ideal case, where only hydrodynamic forces are present and the Reynolds num- 
ber is exactly zero, the motion of the particles is deterministic and reversible due to 
the linearity of the governing flow equations, thereby implying that, upon reversing the 
direction of flow, the particles should retrace their trajectories. However, in any physical 
experiment, there exists inherent slight irreversible effects at the microscopic level, such 
as residual Brownian motion and surface-roughness effects, which, as is usually the case 
in dynamical systems, can have a strong impact on macroscopic measurable quantities. 
As an example, let us consider again the loss of fore-aft symmetry in sheared suspensions. 

It can be shown, based upon the reversibility of Stokes flow, that the pair distribution 
function of sheared suspension of perfect spheres should have fore-aft symmetry. On 



the other hand, since the original work by Gadala-Maria & Acrivos (1980), there exists 
strong experimental evidence that, even at vanishingly small Reynolds numbers and 
Brownian force effects, sheared suspensions develop an anisotropic structure resulting 



in the loss of fore-aft symmetry (Parisi & Gadala-Maria 1987; Rampall et al. 1997). 
This broken symmetry has been attributed to surface-roughness effects. As mentioned in 



2.1 
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the minimum separation between two colliding spheres in a shear flow can be less 



than 10 4 of their radius, so that even a small surface roughness can have an important 



influence in this case (da Cunha & Hinch 1996). Thus, small irreversible effects related to 
surface roughness, which is present at microscopic scales, have a measurable impact on 
the macroscopic structure of the suspension and correspondingly on related macroscopic 
quantities such as the normal stress difference in sheared suspensions. 

Thus far we have discussed the microscopic origin of irreversibility and its manifesta- 
tion in macroscopic quantities which provides a microscopic basis for the existence of an 
intrinsically irreversible macroscopic description as given by the diffusion equation. But, 
the basic assumption underlying the derivation of such an equation and, in particular, 
the validity of a statistical description of the system based on the randomness of the 
microscopic motion of the particles, is in need of further discussion. Recall that, in the 
calculation of the diffusivity in very dilute sheared suspensions, it is generally assumed 
that successive collisions between spheres are statistically independent, sometimes re- 



ferred to as molecular chaos (Acrivos et al. 1992; Wang et al. 1996, 1998; da Cunha & 



Hinch 1996 ). Furthermore, it has been suggested that a close connection exists between 
molecular chaos and dynamical chaos, according to which stochastic- like behaviour is 



possible even for deterministic mechanical systems (e. g. Gaspard (1998 ), p. 225; Dorf- 



man (1998 ).) In fact, even low-dimensional deterministic dynamical systems are known 



to give rise to diffusive transport (usually named deterministic diffusion (Gaspard 1998, 
p. 293)). Chaotic systems have the property that any small perturbation in the state 
of the system will grow exponentially in time, to the point where the evolution can no 
longer be accurately predicted. But, although the hypothesis of chaotic motion as the 
basic mechanism responsible for this loss of memory and for that matter, as the cause of 
the phenomenon of shear induced diffusion, has been suggested in the context of sheared 
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suspensions by Marchioro & Acrivos (2001) and has been supported by their numerical 



simulations showing irreversible behaviour upon reversal in the direction of the shear, 
to-date the presence of chaos in a sheared suspension has not been examined. We shall 
therefore numerically investigate the presence of chaotic motion in sheared suspensions 
by evaluating the largest Lyapunov exponent which is a standard measure of chaoticity 
( Schuster 1989j P- 24). 



It is worth mentioning here that we do not wish to imply that diffusive motion cannot 
possibly operate in the purely hydrodynamic case, and that diffusion owes its existence 
to the unavoidable presence of small microscopic irreversible forces. In fact, the presence 
of chaotic motion, which we shall presently demonstrate, strongly suggests that, even 
in the absence of such effects, diffusive behaviour should still arise due to the loss of 
correlation in the particle motions. 

3.1. Largest Lyapunov exponent 

In a sheared suspension at zero Reynolds number, the state of the system is fully de- 
termined by the coordinates of all the particles and, therefore, a sheared suspension 
consisting of N particles can be considered as a dynamical system in a 3-ZV-dimensional 
phase space. A point of this phase space T is given by 37V particle coordinates and the 
Lyapunov exponents measure the average rate of separation of two initially neighboring 
trajectories in phase space. Thus, given two states of the system separated a small dis- 
tance d(0) in phase space, the largest Lyapunov exponent (LLE) for that state is formally 
defined as 

'1 d(t)' 



A = lim lim 

t->oo d(0)^0 



(3.1) 



t d(0)_ 

with d(t) being the separation distance in phase space at time t, which controls the 
exponential divergence of initially close trajectories. 
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Relaxation to steady state Separation in phase-space 

Figure 4. Schematic representation of the generation of initially close trajectories in a stationary 

state. Starting from a random configuration of hard-spheres I\ we follow the evolution of the 
system towards steady state. Once in steady state, we generate a slightly perturbed state T p by 
adding a small random displacement to each particle. The initial separation in phase space is 
d(0). Then we follow the evolution of the two systems by computing the distance at each time 
d(t). 

In order to obtain a numerical estimate of the LLE for the steady state of a sheared 
suspension, we compute the ensemble averaged LLE used in ergodic systems ( |Gaspard 



1998, p. 144, 247) by generating a number N c of initial configurations T s (0) in a steady 



state flow, starting from a random configuration ITj of spheres and evolving the system 
for a sufficiently long time (typically for strains t ~ 50). For each of these configurations 
we then introduce a slightly perturbed state r p (0), by adding a random displacement d 

to r s (0)|(see figure @). 

_f In order to generate a random perturbation in phase space d, with constant magnitude 
\d\ = d(0), we initially construct a 3N-dimensional vector £, where each component £i is a 
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Figure 5. Separation of initially close trajectories in phase space d(t) for different initial states 
of the unperturbed system. The results correspond to simulations with N = 64, Fo = 1.0, 
r c = 10" 4 , and <f> = 0.35. 

Following the evolution of these two independent systems we can compute a distance 
in phase space given by, 

d{t) = \\T s {t) -T p (t)\\ = (j2 -x*(t)]^ (3.2) 

where the indices s and p refer to the unperturbed and perturbed initial states of the two 

systems! (see figure ||). To be sure, for a single realization, d(t) will depend on both initial 

states r s (0) and r p (0), but, after some transient behaviour, we should certainly expect 

random variable with uniform distribution in the range —1 < < 1. Then, we renormalize 

this vector to obtain the desired initial magnitude of the perturbation, d — f (d(0)/\£\). Thus, 

a random displacement d in phase space corresponds, on average, to a random perturbation, of 

order d(0)/y/W, to the position of each particle. 

| Let us mention that, even though the afhne shearing motion of the spheres is not removed 
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an exponential separation with large fluctuations due to the details of the dynamics of 

the system. This is illustrated in figure ^|. In order to compute the LLE accurately, we 

therefore first average over N c different initial conditions in phase space, then take the 

logarithm of the mean exponential separation, 

A r (t) = In «d(t)) r ) = In ^ £M*)}*) ( 3 - 3 ) 

and identify Ap (the largest Lyapunov exponent) as the slope of Ap{t) in the region of 
its linear growth |. 

Analogously, we can define a distance in the space of the transverse velocity components 
of all the particles, 



A„(t) = ln 



1 N ° \ ( N 



fe=i 



(3.4) 



However, as d(t) grows exponentially, so does a generic projection in phase space (or a 
distance measured with any particular metric) Moreover, since the separation is domi- 
nated by the exponential growth with the LLE, and the velocity is the derivative of the 
and its contribution dominates the initial growth of the separation distance in phase space, the 

exponential growth should eventually become dominant at long times. 

| The method used in this work to compute the largest Lyapunov exponent differs from 



that described by Benettin et al. (1976) and used frequently. Instead of averaging over different 



realizations (ensemble average), Benettin et al. (1976) reset the distance between trajectories to 
the initial value after fixed intervals of time r and compute the Lyapunov exponent from the 
average distance reached before the rescaling procedure (time average). However, both methods 
should give the same value of A if the system is ergodic, because in that case, the time-average of 
any dynamical quantity is equal to its ensemble average over a large number of realizations N c 



Eckmann & Ruelle (1985). The main advantage of our method is that it allows us to compute the 



evolution of all the initial configurations simultaneously, which reduces the computational time 
enormously, as compared to a single very long simulation, when the computations are performed 
using a cluster of processors. 
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Figure 6. Separation of initially close trajectories in phase space, (a) Coordinate phase space, 
(b) Transverse velocity space. The results corresponds to a sheared suspension of N = 64 parti- 
cles, volume fraction cj> = 0.35, interparticle force Fo = 1.0, characteristic range r c = 10 -4 , and 
iV c = 100. The solid line corresponds to a linear fit with (a) A = 0.55 ± 0.04, (b) A = 0.56 ± 0.04. 



The dashed line in (b) is the asymptotic value of A„ as t — *• oo given by (3.5) 
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position of the particles, one would expect the Lyapunov exponent X v , as obtained from 

the slope of A 1 ,(t), to be the same as that measured from Ar(i), i.e. A = Ar = A„. 

In figure |^ we show the time evolution of Ap (t) and A„ (t) for a sheared suspension with 
N = 64 particles, volume fraction <f> — 0.35, interparticle force F — 1.0 and characteristic 
range r c = 10~ 4 . The distance d(t) was averaged over iV c =100 different initial conditions. 
The initial distance in phase space T is d(0) = 10 -4 , which corresponds to a random 
displacement ~ 10 -5 added to the position of each particle. Simulations using a smaller 
initial displacement give a similar behaviour with a variation in the measured Lyapunov 
exponent less than 5%. We also performed simulations with r c = 10 -3 as well as with 
Fq = 0.1, but, in all cases, the variations between the LLE's were within 10%. However, 
larger changes should be expected in the value of the LLE for larger variations in the range 
or the strength of the interparticle force. Finally, the Lyapunov exponent calculated from 
Ar(t) and from A v (t) is, as expected, the same, within the error in its determination. 

In figure ^ three different regimes can be directly observed, a short initial transient be- 
haviour, a large linear growth corresponding to the exponential divergence in phase space 
and a deviation from the linear growth at long times, corresponding to an asymptotic 
saturation in velocity space. 

There are two related effects that might cause the short-time transient behaviour. On 
one hand, many-body nonlinear effects only come into play after at least one collision 
has taken place, and we can therefore relate the transient time to a characteristic time 
between collisions. During this initial time the distance grows mainly due to the shear 
flow imposed to the particles. However, the definition of such a characteristic time is 
not clear, given the fact that the hydrodynamic forces are long ranged. The experiments 
of Brccdveld et al. (2001 6| ) found a similar transition when measuring the diffusivity of 
the particles. They suggested that this transition might correspond to the deformation, 
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Figure 7. Dependence of the Lyapunov exponent on the volume fraction. The results correspond 
to numerical simulations for N = 64, N c — 100, Fq = 1.0 and r c — 10~ 4 . Ar is computed from 
the separation growth in coordinates phase space, and X v using the transverse velocity space. 

due to the shear flow, of the initial cloud of particles surrounding a test sphere. In our 
simulations we found that this transient time decreases with increasing concentration, 
which is consistent with the suggestion referred to above. 

In the simulations it is observed that the distance in the transverse velocity space 
saturates. This saturation takes place because the definition of the distance in the trans- 
verse velocity space does not include the affine motion of the particles, and therefore, two 
randomly chosen systems in steady state will have a finite mean distance. In fact, the sat- 
uration distance can be estimated by assuming that the perturbed system loses memory 
of its initial unperturbed state. In this case, if we further assume that A£° becomes 



A° 



In 



N v! 



N 



In [27V (a 2 vy +a 2 vz )} 



(3.5) 
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where <r^ is the variance in the i-th component of the velocity, the saturation constant 

can be computed in the steady state, independently of the evolution. For the simulation 

shown in figure^, the variances in the velocity components are o~% y — 0.20 and a^ z = 0.07 

respectively, giving a saturation value = 1.77 which is in excellent agreement with 

the observed saturation distance, as shown by the dashed line in figure ||. On the other 

hand, as the velocity-space distance saturates, the system loses its memory of the initial 

state so that the separation in coordinate-phase deviates from the exponential growth, as 

can be observed in figure |[ In fact, the loss of memory of the initial state of the system 

leads to an asymptotic diffusive behaviour in phase space. 

3.2. Dependence of the Lyapunov exponent on the concentration 

We investigated the dependence of the largest Lyapunov exponent on the volume fraction 
of the suspension. As the concentration increases, the frequency of collisions also increases 
and we may expect an enlarged sensitivity to the initial conditions. In figure ^ we show 
the numerical values of the LLE as a function of <j>. It can be seen that the stochasticity of 
the system increases with concentration, as measured by A. An almost linear dependence 
can be also observed, with the surprisingly fact that the LLE appears to remain finite as 
0^0. We remark parenthetically that, in kinetic theory the Lyapunov exponent can be 



shown to be a linear function of the frequency of collisions ( Gaspard (1998 ), p. 8; Dorfman 



(1998 )), and that in a sheared suspension, the frequency of collisions is roughly v ~ <f>. 
However, due to the long-range interactions between the particles, it is not clear that 
the same relation between the LLE and v will continue to apply in the case of sheared 
suspensions; this then is the main difference between the present situation and kinetic 
theory, which is based on the existence of short-range interactions between the particles. 
In fact, the presence of these long-range hydrodynamic forces might be responsible for 
the curious behavior, seen in figure |7|, of A at small volume fractions. Unfortunately, the 
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observed initial transient time increases as the concentration becomes smaller, thereby 

limiting the range of concentrations which we were able to investigate to (f> > 0.01. Thus, 

the functional dependence of the LLE on <f> in the dilute limit remains an open problem. 



4. Transverse velocity autocorrelation function 



In the previous section we discussed the global loss of memory of a system when 
its initial state is slightly perturbed. At the level of individual particles this asymptotic 
stochastic motion is also present, and could be studied by examining whether the velocity 
of a particle decorrelates from its initial value. The average correlations between the 
velocity of a particle at two different times are measured via the velocity autocorrelation 
function. Here, we consider K a (t), where a refers to the transverse velocity component 
being measured 



K a (t) = (v a (0)v a (t))/a 2 Va 



c fe=i l i=l 
and its time integral r Va 



1 Nc { 1 N 



c fc=i I, i=i 



>v a JO 



dt(v a (0)v a (t)) = / dtK a (t) 



,(4.1) 



(4.2) 



The importance of these functions in computer simulations is well known (Allen & 



Tildcsley 1987, p. 58) and, in particular, the velocity autocorrelation function is related 
to the diffusion coefficient through, 



(v a (0)v a (t))dt = dl T Va 



(4.3) 



This expression has been used in the context of suspensions by Nicolai et al. (1995) 



to determine the diffusivity of sedimenting non-Brownian spheres, and by Marchioro & 



Acrivos (2001) to measure the shear-induced self-diffusivities. 



In figure || we show the computed velocity autocorrelation functions for both transverse 
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velocity components, K y (t) and K z (t), for volume fractions ranging between 0.05 < (j> < 

0.45, where y and z are, respectively, along and normal to the plane of shear. In all cases, 

the velocities becomes uncorrelated at long times, around t = 10 for the lowest volume 

fraction, and from then on oscillations around zero could be attributed to statistical noise. 

Thus, for t beyond about 10, the displacement of an individual particle is independent of 

its displacement at previous times and hence should be describable in terms of a random 

walk. 

It can also be observed that, except at high concentrations (0 w 0.45), the autocorre- 
lation function becomes negative over a range of several time units. Let us mention, that 
a similar behaviour for the velocity autocorrelation function was found in molecular dy- 
namics simulations of simple liquids ( Alder fc Wainwright 1958| ; Rahman 1964), and was 



explained in terms of the backscattering by neighboring particles at high densities ( Aider 



& Wainwright 1967, 1970). In fact, the negative correlation in velocity suggests that, on 
average, individual particles reverse their velocity, typically after a time interval of order 
1 according to our simulations. This result can be understood in terms of the motion 
of colliding particles in linear shear flow. Specifically, as was already mentioned in the 
introduction, when two isolated spheres collide, the net displacement in both transverse 
directions is zero, as a consequence of the reversibility of the creeping flow equations and 
the symmetry of the problem. However, the instantaneous deviation in the velocity of 
the spheres during the encounter is not zero and clearly anti-correlated, since the net 
displacement should integrate to zero. Thus, encounters between two isolated spheres 
undergoing linear shear flow give rise to a negative correlation in the velocity fluctua- 
tions. On the other hand, encounters involving three or more particles will, in general, 
not be symmetric and the particles will experience a net displacement from their origi- 
nal streamlines ( Wang et al. 1996| ). Therefore, in general, the interaction between more 
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Figure 9. Average number of particles with their transverse velocities reversed from their initial 
direction. Different curves corresponds to different components of the velocity as indicated. The 
simulations are for N = 64, N c = 100, F = 1.0, r c = 10~ 4 and <j) = 0.35. 

than two particles does not necessarily contribute to a negative autocorrelation in the 
transverse velocity. 

The previous discussion not only explains the observed negative correlation but also 
the fact that it becomes more pronounced with decreasing volume fraction and that, 
as <j> — > 0, it seems to converge to an asymptotic function dominated by two-particle 
encounters. 

At high concentrations we see that both velocity correlation functions decay rapidly 
to zero and show very little structure. This effect might be attributed to 'screening': i. e. 
each of the particles alters the ambient velocity field, and when these perturbations have 
a high spatial density and a rapid time dependence, their effect is to produce a strongly 
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fluctuating background flow which tends to decorrelatc the velocity of any particle from 

its earlier value. An analogous effect is present in porous media flows, which could be 

thought of as an extreme case of a suspension so concentrated as to be immobile. Here, a 

static force perturbation in a fluid-saturated porous medium decays exponentially with 



distance (as seen from the Brinkman equation ( |Brinkman 1947| ; |Dmiofsky fc Brady 1987| ), 
for example), in contrast to the power-law decay in more dilute mobile suspensions. 

In figure [)] we present the time evolution of the average fraction of particles having 
their initial velocities reversed, 

V a (0)v a (t) 



n(t) 



1 - 



(4.4) 



K(o)v a (t)\ 

It can be seen that, on average, more than half the particles have both transverse velocity 
components reversed at intermediate times (t ~ 1) for the indicated values of Fq, r c and 
<f). These results are consistent with our previous discussion regarding the origin of the 
negative autocorrelation. 

Finally, let us discuss some of the important consequences of the negative correlations 
in the transverse velocity and their origin. In the first place, it is clear that the integral 
time scale given by equation (12) is different from the autocorrelation timqjj This is 



illustrated in figure ttOj, where we depict the values for r Vy and t Vz , obtained using (4.2), 
as a function of the volume fraction. It is clear that the integral time scale measured 
from (fh^) fails to capture the actual time scale underlying the loss of correlation in 
the particle velocities in that it increases with concentration, contrary to the fact that 
the actual time scale for the correlations becomes shorter as the concentration and the 
rate of collisions between the particles increases. Also, given the limited accuracy in the 



t A correlation time denned through equation (4.2) implicitly assumes an exponentially de- 



caying autocorrelation function, as is the case for any Gauss-Markov process ( ra,n Kampen 1987 , 
p. 81). 
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Figure 10. r Vy and t Vz given by equation ( |4.2| ), as a function of the volume fraction. The 
simulations are for N = 64, N c = 100, F = 1.0, and r c = 10" 4 



computation of the integral (4.2), the values of r Vy and t Vz shown in figure |10| are subject 
to large relative errors when <\> < 0.10. 

The same difficulty is encountered when the diffusion coefficient is evaluated by inte- 
grating the velocity autocorrelation function. As previously discussed, the motion gen- 



erated by binary collisions is not diffusive and thus, in view of equation 4.3, the integral 
of the autocorrelation function should vanish in the limit <j) — ► 0. This is born out by the 
results shown in figure |l0| where it can be seen that the integral of the autocorrelation 
function approaches zero with decreasing volume fraction. Thus, the leading contribution 
to the diffusivity comes from a small portion of the autocorrelation function the role of 
which becomes increasingly more crucial as the volume fraction decreases. 
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5. Velocity fluctuations 

In the previous sections we discussed the loss of correlation in the particle velocities 
and how the integral time scale is related to the diffusivity. On the other hand, as 



is clear from equation (4.3), the magnitude of the velocity fluctuations also plays an 
important role in determining the particle diffusivity. Understanding velocity fluctuations 
in the presence of long-range hydrodynamic interactions is a long standing problem, 



and has been recurrently studied in the context of sedimcnting suspensions (Brenner 



fc Mucha 2001 ). In this section we shall present numerical results concerning not only 
the magnitude of the velocity fluctuations given by a v but also the whole probability 
distribution function (pdf ) of the velocity fluctuations, for different volume fractions. 

Numerical simulations clearly suffer from two limitations, which become of particular 
importance when computing velocity fluctuations: the small number of particles in the 



simulations and the necessary approximations to the hydrodynamic forces (Brady & 



Bossis 1988). However, we shall show that the functional form of the pdf's undergoes a 



clear transition, from an exponential to a Gaussian distribution as the volume fraction 
is increased, which is not sensitive to the number of particles. Moreover, this feature 
could be measured experimentally thereby providing a useful test for the approximations 
involved in the hydrodynamic interactions which have been used in constructing the 
Stokesian dynamics code. 

In figure [□] we show the pdf of the velocity fluctuations in the direction of the shear 
P(v y ), obtained at two different volume fractions, <fi = 0.05 and <f> — 0.35, and using 
different number of particles in the simulations. It can be seen that, as previously stated, 
the distribution of velocities is insensitive to the number of particles, at least as far as 
its functional form is concerned. 

In figure O we compare the pdf's at low and high volume fractions in a log-linear 
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Figure 1 1 . Probability density function of the velocity fluctuations in the direction of the shear 
P{v y ) for two different volume fractions 4> = 0.05 and <f> = 0.35. Different curves correspond 
to different number of particles in the numerical simulations. The simulations are for N = 64, 
N c = fOO, F = 1.0 and r c = f0~ 4 . 
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Figure 12. Probability density function of the velocity fluctuations in both transverse directions 
(a) P(v y ) and (b) P(v z ), for two different volume fractions (f> = 0.05 and (f> = 0.35. See §|B| for a 
discussion of the different distributions used to fit the numerical data (solid and dashed lines). 
The simulations are for N = 64, iV c = 100, F = 1.0 and r c = 10" 4 . 
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plot. Clearly, both P(v y ) and P(v z ) become not only broader as the concentration 



increases but change noticeably in their shape. In figure |12ja we show that, at low 
volume fractions {<p — 0.05), P(v y ) is well described by an exponential distribution, 
P(v y ) oc exp(— 11.4 \v y \) but that, at a large volume fraction, (f> = 0.35, a Gaussian distri- 
bution accurately fits the numerical results, P(v y ) oc exp(— 2.41 Vy). A similar behaviour 
is observed in P(v z ), as shown in figure p~2| b in that, at low volume fractions, the pdf is 
well approximated by an exponential distribution P(v z ) oc exp(— 26.3 \v z \), but at large 
concentrations neither a Gaussian nor an exponential distribution accurately fits the 
numerical results over the whole range of velocities. However, we show that, for small 
velocity fluctuations, the distribution is Gaussian (dashed line; P{v z ) oc exp(— 7.76 v%)), 
whereas large fluctuations are better described by an exponential distribution (dotted 
line; P{v z ) oc exp(— 7.57 |i> z |)). A similar transition, but in the opposite direction, was 
observed in fludized suspensions of non-colloidal particles, where the velocity fluctuation 
distributions vary from Gaussian to exponential as the particle concentration increases 



( |Rouyer et al. 1999| , gOOOj ). 

We believe that this transition in the pdf of the velocity fluctuations towards a Gaus- 
sian distribution, is due to the predominant role played by lubrication forces at the 
high concentrations. Recall that, given their short-range character, lubrication forces are 
essentially two-body interactions and are accounted for in a pairwise additive way in 
the simulation method. On the other hand, again in the simulation method, many-body 
long-range hydrodynamic interactions are accounted for by means of the far-field approx- 
imation^ Therefore, the random addition of lubrication forces, generated from spheres 

f A detailed discussion of the approximations involved in the Stokesian dynamics method 



can be found elsewhere (Durlofsky et al. 1987; Brady & Bossis 1988; Brady et al. 1988; 



Ichiki 



Brady 2001). 
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in close proximity to one another, is the dominant contribution at high concentrations, 

and would be expected to yield a Gaussian distribution. Let us also note that the hydro- 
dynamic forces inhibit the relative motion of the particles as their approach one another 
and therefore would prevent large fluctuations in the velocity from occurring, consistent 
with the range in v z where a Gaussian distribution properly describes the corresponding 
pdf's. On the contrary, at low concentrations, lubrication forces becomes negligible and 
therefore long-range many-body interactions determine the fluctuations in velocity. In 
this case, our results show an exponential distribution. It is clear then, that a compari- 
son with experimental results would provide another test of the numerical method and 
of the approximations that are involved in the simulations. 

An alternative view of the distinction between Gaussian and exponential velocity distri- 
butions may be offered by analogy with turbulence. Specifically, in boundary layer flows, 
one observes that the pdf of vorticity is Gaussian at low Reynolds numbers, but develops 
exponential tails at the very high Reynolds numbers characteristic of atmospheric flows 



(Fan 1991). Similarly (perhaps), the pdf of temperature fluctuations in thermal convec- 
tion exhibits a transition from Gaussian to exponential as the Rayleigh number increases 
( [5ano et al. 198S ; Wu 1991 ). In both cases, the transition is gradual rather than abrupt, 
and the increase of Reynolds or Rayleigh number is accompanied by the development 



of organized large-scale coherent structures in the flow ( Frisch 1995 , p. 100). Since one 
normally associates Gaussian pdf's with the addition of uncorrelated random variables, 
it is natural to relate the transition to the more slowly-decaying exponential pdf's to 
the appearance of correlated long-range structures. In the case of suspensions, inverse 
concentration is analogous to the Reynolds or Rayleigh number in the sense that a large 
value of the parameter is associated with long-ranged correlations, as we have seen in 
the discussion of screening in §^|. A common feature in these three situations is that, 



Dynamics of sheared suspensions 31 
at low values of the appropriate control parameter, the motion is relatively uncorrelated 

and the fluctuations are Gaussian while, at high values, the motion is organized and the 

fluctuations are more persistent and exponentially distributed. 



6. Shear-induced self-diffusion at low concentrations 

In section §|| we discussed how stochastic (diffusive- like) transport arises in the deter- 
ministic dynamics of sheared suspensions, due to the chaotic evolution in phase space 
and, in we showed the presence of a stochastic motion at the individual level of a 
single sphere by studying the loss of correlation in the transverse particle velocities which 
leads to a random diffusive motion of single particles. Finally, in figure we show that 
the mean square displacement of a single particle becomes diffusive after a time approx- 
imately equal to the characteristic time after which the particle's velocity correlation is 
lost. 

We also discussed the relation between the negative correlation in the transverse ve- 
locity fluctuations and the hydrodynamic interaction between a pair of non-Brownian 
spheres undergoing shear. We also remarked that, in order to get a diffusive motion from 
purely hydrodynamic interactions, collisions between at least three particles are necessary 
while, on the contrary, in the presence of non-hydrodynamic forces such as the interpar- 



ticle repulsive force introduced in § 2.1, binary collisions alone lead to diffusive motion. 
This fact gives rise to different scaling relations for the diffusivity depending on the rela- 
tive importance of the interparticle force and the hydrodynamic forces. Let the diffusion 
coefficient in the pure hydrodynamic limit be denoted by D^, and the contribution to the 
diffusivity in the presence of non-hydrodynamic forces by D n _h. As already mentioned 
in §[l], since the hydrodynamic contribution arises from collisions between three or 
more particles, and the rate at which a given sphere interacts with two other spheres is 
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Time 

Figure 13. Mean square displacement ol a single particle, averaged over all the N = 64 particles 
in the suspension and over N c — 100 diHerent realizations, as a function of time for different 
volume fractions. The magnitude of the interparticle force is Fo = 1 and the characteristic range 
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Figure 14. Dimensionless diffusion coefficient in the direction of the shear D y as a function of 



the volume fraction. The solid line correspond to the theoretical result calculated by Wang et al 



(1996 ), in the absence of non-hydrodynamic forces. The simulations are for N = 64, N c = 100, 
Fo = 1.0 and r c = 10~ 4 . 

proportional to 7 cf> 2 in the limit of low volume fractions, should scale as 7 (j> 2 a 2 as 
(f> — ► while D n _h oc f{Fo, r c ) j(pa 2 , where f(Fo, r c ) should be an increasing function of 
the interparticle force magnitude Fq, and of the range r c . Then, when interparticle forces 
are negligible, the diffusion coefficient D will be dominated by the hydrodynamic contri- 
bution and should scale as j(j> 2 a 2 . On the other hand, at low enough concentrations or 
large interparticle forces, binary collisions become dominant and a linear dependence on 
the volume fraction should be expected, i. e. D ~ D n -h oc j<fia 2 . 

Let us first investigate the case when the interparticle force is negligible. In figure 14 
we present, in a log-log plot, D y the diffusion coefficient in the direction of the shear 
rendered dimensionless by 7a 2 , as a function of the volume fraction where, in all cases, 
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Figure 15. Dimensionless diffusion coefficient D z as a function of the volume fraction. The 
solid line corresponds to a fit D z — 0.07 4> 2 . The dashed line corresponds to the result obtained 



by [Wang et al. (1996| ), D z = 0.005 cp 2 . The simulations are for N = 64, N c = 100, F = 1.0 and 
r c = 10~ 4 . 

the diffusion coefficient is obtained from the slope of the mean square displacement in 
the region of its linear growth, as shown in figure [l3| In the numerical simulations the 
characteristic range of the interparticle force was set to r c = 10~ 4 and Fo = 1.0. It can 
be observed that at low volume fractions, a quadratic dependence of D y on cj> is obtained, 
as expected when the non-hydrodynamic force is negligible. We also compare in figure 



|l4j the numerical results with the theoretical values obtained by Wang et al. (1996), 
who evaluated the diffusion coefficient by computing the mean square displacement of 
a test sphere over all possible encounters with two other spheres in the absence of the 
non-hydrodynamic force and found that D y = 0.11 </) 2 , in excellent agreement with our 
fitted value D y = (0.11 ± 0.02) <f> 2 . However, it should be kept in mind that the diffusion 
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coefficients reported here are from simulations using only 64 particles, a number which 

may not be large enough for an accurate comparison with other vales of D obtained 

theoretically or experimentally. 



In figure 15 we present the dependence of D z on cf>. In this case, a quadratic dependence 



is also found. However, the theoretical values calculated by [Wang et al. (1996| ) are much 
smaller than our numerical results in that a fit to the numerical results gives D z — (0.07± 



0.007) 2 , whereas Wang et al. [1996 ) found D z = 0.005 <j) 2 . However, let us note that the 
anisotropy found in our simulations for the self-diffusion coefficient, D y /D z ~ 1.5, is in 



agreement with the experimental results of Phan & Lcighton (1999) recently confirmed 



by greedveld et al. (2001a| ). 

Let us now consider the case in which the interparticle force contributes significantly 
to the diffusion coefficient. To this end, we performed simulations with the characteristic 
range of the interparticle force increased 1000 times relative to its previous value r c = 
10~ 4 . The strength of the force was kept constant Fq = 1.0. In figure [l6] we present the 
new numerical results for the diffusion coefficient in the direction of shear, as a function 
of the volume fraction. For comparison we also show the previous results. It is clear that 
the diffusivity deviates from the values found previously and becomes larger than before 
at very low concentrations. As <f> — > 0, a linear dependence on <fi is found to accurately 
fit the numerical data, giving D y = (7.6 ± 1.3) x 10~ 3 <fi. 

In figure [Tt] the diffusion coefficient in the vorticity plane is compared for the same two 
values of the characteristic range, r c = 10~ 4 and r c = 10 , as a function of cf>. In this 
case, the interparticle force is not strong enough and the concentration not low enough 
for the linear regime to be observed. 

It is important to note that the largest range of the interparticle force r c = 0.1 is 
clearly too large to be considered as representing residual Brownian forces or electro- 
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Figure 16. Dimensionless diffusion coefficient in the plane ol shear D y as a function of the 
volume fraction cj>. Open circles corresponds to short-ranged interparticle forces r c = 10~ 4 , 
while solid circles correspond to a longer range r c = 10" 1 . The solid line is a best fit with a 
linear dependence on <f>: D y — (7.6 ±1.3) x 10 -3 cj>. The simulations are for N — 64, N c = 100, 
and F = 1.0. 

static repulsion ( Foss fc Brady 2000| ). However, this force might resemble the effect of 



particle roughness. By means of numerical simulations, da Cunha & Hinch (1996) and 



recently Zarraga & Leighton (2001), showed that particle roughness lead to a diffusivity 
proportional to the volume fraction. These authors modelled particle roughness by a nor- 
mal force which prevents the particles from coming closer than a minimum dimensionless 
separation 2 + e between the centres of the two spheres. As a crude approximation, we 
might estimate the magnitude of the roughness represented by our interparticle force, as 



its characteristic range r c , i.e. e ~ r c ~ 0.1 in our simulations. For e ^ 0.1 da Cunha fc 



Hinch (1996 ) and Zarraga fc Leighton (2001 ) found that D y /(f) ~ 0.02 , which is larger 



Dynamics of sheared suspensions 



37 




in 



10 



4> 



10 



Figure 17. Dimensionless diffusion coefficient perpendicular to the plane of shear D z as a 
function of the volume fraction <j>. Open and solid circles corresponds to different characteristic 
range of the interparticle force, r c = 10~ 4 and r c = 10 _1 respectively. The solid line is the best 
fit with a quadratic dependence on cj>, D z = (0.07 ± 0.007) 4> 2 . The simulations are for N — 64, 
N c = 100, and F = 1.0. 

than our fitted value D y /4> ~ (0.0076 ± 0.0013) . However, as shown in figure |l8|, in 
our simulations, pairs of particles are allowed to come closer than r c .This is because, 



at r = r c the interparticle force is not infinite as in the analysis by da Cunha & Hindi 



(1996) and p^arraga fc Leighton (200l] ). On the other hand, if we consider the observed 



minimum distance between spheres as the effective roughness represented by our repul- 
sive force then, e ~ 0.01 and from da Cunha fc Hinch (1996] ) D y /4> ~ 0.006 , which is 
now in fairly good agreement with our fitted value. A linear behaviour D y oc 4> was also 
found in experiments by Bicmfohr et al. (1993| ) and Zarraga fc Leighton (1999 ). However, 
the diffusion coefficient reported in these experiments was substantially larger than that 
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Figure 18. Pair distribution function g(r). The characteristic range is r c — 0.1. N = 64, 

4> = 0.10, F = 1.0, and iV c = 100. 

obtained by means of numerical simulations {D y /<j> ~ 0.03 as compared to D y /<p ~ 0.005 
from the numerical simulations for e ~ 10 -3 .). 

7. Summary 

The complex dynamics of homogeneous sheared suspensions of monodisperse, neutrally 
buoyant, non-Brownian spheres at zero Reynolds number was investigated by means 
of numerical simulations using Stokesian dynamics. Starting from a large number of 
independent initial configurations (N c = 100), the evolution of typically N — 64 spheres 
undergoing simple linear shear was simulated during a time t ~ 100, which was sufficiently 
long to allow us to study the dynamics of the system in steady state. In addition to 
the hydrodynamic interactions between spheres, the simulations included a short-ranged 
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repulsive interparticle force that qualitatively models the effects of surface roughness and 

Brownian forces both of which play an important role when neighboring spheres nearly 

touch one another. 

We began by recalling some of the well known effects of the non-hydrodynamic in- 
teraction on the microscopic structure of sheared suspensions ( |Bossis fc Brady 1984 ; 



Rampall et al. 1997|) . We showed that the location of the first peak of the pair distribu- 



tion function strongly depends on the range of the interparticle force r c in that, as r c is 
increased, the minimum separation between neighboring spheres increases significantly. 
We also discussed how the angular distribution of pairs depends on this force, showing 
that, as the range of the interparticle force is increased, the pair distribution function 
becomes increasingly more asymmetric, with fewer pairs being oriented downstream than 
upstream. This asymmetry implies the loss of time reversibility and therefore identifies 
the interparticle force as a microscopic origin of irreversibility. 

The dynamics of sheared suspensions was shown to be chaotic over the whole range 
of volume fractions investigated, 0.01 < <f> < 0.35, with the largest Lyapunov exponent 
increasing linearly with cf>. The existence of chaos in the dynamics of the suspensions 
provides an explanation for the loss of correlation in the particle motions leading to 
diffusive behaviour at long times. In fact, we showed that the system loses memory of 
its initial state exponentially, and that the asymptotic separation distance between two 
initially close trajectories in the velocity-space of the two transverse velocity components 
can be accurately estimated assuming that the two systems are completely uncorrelated. 

At the level of individual smooth spheres, we also demonstrate the loss of memory 
of their instantaneous velocity fluctuations via the transverse velocity autocorrelation 
function. For high concentration, ip ~ 0.40, a screening-type mechanism is observed and 
the correlations in particle velocities are lost after a short period of time of 0(1). On 
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the other hand, at the lower concentrations, the transverse particle motions remain, on 

average, correlated for a longer period of time. We found furthermore that, for both the 
transverse components of the velocity, the autocorrelation function becomes negative at 
intermediate times (t ~ 2). We proposed an explanation for this effect based in the dy- 
namics of two isolated spheres undergoing simple shear, according to which, since the 
purely hydrodynamic interaction between two spheres does not lead to any net lateral dis- 
placement, it is a source of negative correlation in the velocity fluctuations during binary 
collisions. This explanation is consistent with the fact that the region within which the 
velocity autocorrelation function is negative continues to enlarge with decreasing <f) with 
the whole autocorrelation function appearing to converge to an asymptotic distribution 
dominated by binary interactions between spheres. We mentioned that the integral of this 
asymptotic distribution must vanish on account of its being proportional to the net dis- 
placement experienced by a pair of spheres, and showed that, in fact, the numerical value 
of the integral of the autocorrelation function steadily decreases as (f> — ► 0. An important 
consequence is that the leading contribution to the diffusivity of the particles as — > 
comes from an asymptotically negligible contribution to the autocorrelation function. 
Therefore, an estimate of the diffusion coefficient, based on the velocity autocorrelation 
function, will be highly inaccurate at least for low volume fractions. 

At this microscopic level we also computed the velocity probability distribution func- 
tion in both lateral directions and observed a transition from an exponential to a Gaussian 
distribution as the volume fraction is increased. We proposed that this transition is due 
to the dominant role, at high concentrations, of the lubrication forces, which are essen- 
tially two-body random interactions and therefore would be expected to give rise to a 
Gaussian distribution. Unfortunately, there are no experimental measurements thus far 
of the probability distribution function of the velocity fluctuations. 
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Finally, we investigated the scaling of the diffusion coefficient D at low concentrations 

as the range of the interparticle force is varied, by setting D equal to one half the slope 

of the mean square displacement in the region of its linear growth. For a very small 

range, r c = 1CP 4 , we showed that the effect of the interparticle force is negligible and D 

was found to scale as 7</> 2 a 2 for 0.03 < <j> < 0.15, both along and normal to the plane 

of shear. This results corresponds to a diffusive motion arising from collisions between 

three spheres simultaneously and, in the case of D y , is in very good agreement with 



the values obtained by Wang et al. (1996 ), who evaluated the diffusion coefficient by a 
completely different procedure, viz. by computing the mean square displacement of a test 
sphere over all possible encounters with two other spheres. On the other hand, for a much 
larger range of the interparticle force, r c = 10 _1 , a linear dependence of the diffusion 
coefficient on <f> is observed as <j) — » 0, which implies that a significant contribution to the 
diffusion coefficient emanated from encounters between two particles. This is consistent 
with the notion that the interparticle force qualitatively mimics the effects of surface 
roughness and other non-hydrodynamic forces, and reproduces the scaling behaviour 
found in the theoretical analysis by |da Cunha fc Hinch (1996 ) and in the experimental 



work by Zarraga fc Leighton (1999| ). However, the numerical simulations still fail to 



reproduce the experimental values of the diffusivity reported by Zarraga & Leighton 



1999 ) by an order of magnitude. 
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